Generalized Dissimilarity Modelling (GDM)

Genomic offset predictions

Published

January 8, 2025

1 Introduction

Resources on Generalized Dissimilarity Modelling (GDM):

From the GDM website: ‘The R package gdm implements Generalized Dissimilarity Modeling (Ferrier et al. 2007) to analyze and map spatial patterns of biodiversity. GDM models biological variation as a function of environment and geography using distance matrices – specifically by relating biological dissimilarity between sites to how much sites differ in their environmental conditions (environmental distance) and how isolated they are from one another (geographical distance). […] GDM also can be used to model other biological levels of organization, notably genetic (Fitzpatrick and Keller 2015) [..] and the approaches for doing so are largely identical to the species-level case with the exception of using a different biological dissimilarity metric depending on the type of response variable.’

2 Formatting data

2.1 Genomic data

2.1.1 Pairwise \(F_{ST}\) matrices

The GDM analysis takes as inputs matrices of pairwise \(F_{ST}\).

To estimate the pairwise \(F_{ST}\), we use the individual-level (i.e. allele counts for each genotype) genomic data with missing data (i.e. no imputation of the missing data), and without minor allele frequencies.

We use the Weir and Cockerham method, with the function pairwise.WCfst of the hierfstat package. We could have used the function gene.dist (with the option WC84), as in Gougherty, Keller, and Fitzpatrick (2021) (and the preprint Gougherty et al. 2020). gene.dist keeps only the lower triangle of the \(F_{ST}\) matrix, while pairwise.WCfst keeps the whole matrix.

Code
# we load the allele counts of each genotype
geno <- read.csv(here("data/DryadRepo/FormattedFilteredGenomicData_AlleleCounts_withoutmaf.csv"), row.names = 1)

geno[geno ==1] <- 12
geno[geno ==2] <- 22
geno[geno ==0] <- 11

geno <- geno %>% 
  t() %>% 
  as.data.frame() %>% 
  dplyr::mutate(pop=substr(row.names(.), 0, 3)) %>% 
  dplyr::select(pop, everything())

geno[1:10,1:10] %>% kable_mydf(boldfirstcolumn = T)
pop snp_2 snp_3 snp_4 snp_5 snp_6 snp_7 snp_8 snp_9 snp_10
ALT10 ALT 11 11 11 11 11 11 11 11 11
ALT2 ALT 11 11 11 11 11 12 11 11 11
ALT3 ALT 11 11 11 11 11 11 12 11 11
ALT4 ALT 11 11 12 11 11 12 11 11 11
ALT5 ALT 11 11 11 11 11 11 11 11 11
ALT6 ALT 12 11 12 11 11 11 11 11 11
ALT8 ALT 11 11 11 11 11 12 11 11 11
ALT9 ALT 11 11 11 11 11 11 11 11 11
ARM1 ARM 11 11 11 11 11 11 11 12 11
ARM10 ARM 11 12 11 11 11 11 11 11 11

We load the sets of candidate and control SNPs.

Code
snp_sets <- readRDS(here("outputs/list_snp_sets.rds"))

We calculate the pariwise \(F_{ST}\) matrices.

Code
fst_matrices <- sapply(snp_sets, function(x){ 
  
geno %>% 
  dplyr::select(pop,all_of(x$set_snps)) %>%
  pairwise.WCfst(diploid=TRUE)
    
}, USE.NAMES = TRUE,simplify=FALSE)

# save it
saveRDS(fst_matrices,file=here("data/GDManalysis/FstMatrices.rds"))

2.1.2 Negative values

Pairwise \(F_{ST}\) matrices contains some negative values.

Code
neg_vals <- sapply(fst_matrices, function(x){ length(x[which(x <0)])}, USE.NAMES = TRUE,simplify=FALSE)

lapply(names(snp_sets), function(x){
  
  tibble('SNP set' = snp_sets[[x]]$set_name,
         'Nb of negative values' = neg_vals[[x]])
  
}) %>% 
  bind_rows() %>% 
  kable_mydf()
SNP set Nb of negative values
All candidate SNPs (380) 6
Common candidate SNPs (69) 34
Candidate SNPs considering pop. struct. correction (221) 8
Control SNPs unmatching allele frequencies (380) 6
Control SNPs matching allele frequencies (380) 8
All SNPs (9817) 2

It generally means there is more variation within than among populations and is likely to result from uneven sample sizes.

2.1.3 Scaling the \(F_{ST}\) matrices

We have to scale the \(F_{ST}\) matrices between 0 and 1 to facilitate model convergence and to enable comparisons the different sets of SNPs, which display different ranges of observed \(F_{ST}\) values.

We standardize the \(F_{ST}\) values in the following way:

\[x_{new} = \frac{x- x_{min}}{x_{max}-x_{min}} \]

This was proposed here: https://www.statisticshowto.com/probability-and-statistics/normal-distributions/normalized-data-normalization/

Code
fst_matrices <- sapply(fst_matrices, function(x){
  
 (x-min(x,na.rm=T))/(max(x,na.rm=T)-min(x,na.rm=T))
  
}, USE.NAMES = TRUE,simplify=FALSE)


# To check that it worked
# sapply(fst_matrices, function(x){range(x,na.rm=T)}, USE.NAMES = TRUE,simplify=FALSE)

2.1.4 Matrix visualization

Code
source(here("scripts/functions/corpmat.R"))
col <- colorRampPalette(c("#BB4444", "#EE9988", "#FFFFFF", "#77AADD", "#4477AA"))

viz_matrix <- function(x,plot_title) {

  p.mat <- corpmat(x)

  p <- corrplot::corrplot(x, method="color", col=col(200),
                   type="upper", order="hclust", 
                   addCoef.col = "black", # Add coefficient of correlation
                   tl.col="black", tl.srt=23, #Text label color and rotation
                   # Combine with significance
                   p.mat = p.mat, sig.level = 0.05, insig = "blank", number.cex =0.8,tl.cex = 0.8,
                   # hide correlation coefficient on the principal diagonal
                   diag=FALSE,
                   title = plot_title,mar=c(0,0,2,0)) # add an upper margin to see the title
  
  return(p)
}
Code
lapply(names(snp_sets), function(x) viz_matrix(x=fst_matrices[[x]], 
                                   plot_title = snp_sets[[x]]$set_name))

2.1.5 Matrix formatting for GDM

For the gdm package: the distance matrix must have as the first column the names of the populations.

Code
fst_matrices <- sapply(fst_matrices, function(x){
  
x <- x %>% 
    as.data.frame() %>% 
    tibble::rownames_to_column("pop") 
  
}, USE.NAMES = TRUE,simplify=FALSE)


saveRDS(fst_matrices,file=here("data/GDManalysis/ScaledFstMatrices.rds"))
Code
fst_matrices <- readRDS(here("data/GDManalysis/ScaledFstMatrices.rds"))

2.2 Climatic data

We do not have to scale the climatic data before the GDM analysis, as scaling of predictors is part of model fitting in GDM, as said in Mokany et al. (2022): ‘As different predictors are measured on different scales (e.g., temperature in degrees, precipitation in mm), they are transformed as part of model fitting, such that the transformed distance between a pair of sites for different predictors can be meaningfully compared and combined.’

Code
# Selected climatic variables
# ===========================
clim_var <- readRDS(here("data/ClimaticData/NamesSelectedVariables.rds"))

extract_climatedt_metadata(var_clim = clim_var) %>% 
  dplyr::select(label,description,unit) %>% 
  set_colnames(str_to_title(colnames(.))) %>% 
  kable_mydf(font_size = 12)
Label Description Unit
bio1 Mean annual temperature Celsius degrees (°C)
bio3 Isothermality (bio2/bio7) (×100) index
bio4 Temperature seasonality (standard deviation ×100) Celsius degrees (°C)
bio12 Annual precipitation Millimeters (mm)
bio15 Precipitation seasonality (coefficient of variation) index
SHM Summer heat moisture index °C / mm
Code
# Loading point estimate climatic data
adj <- "noADJ"  # not adjusted for elevation
ref_period <- "ref_1901_1950" # reference period 1901-1950

clim_past <- readRDS(here(paste0("data/ClimaticData/MaritimePinePops/ClimatePopulationLocationPointEstimates_ReferencePeriods_",adj,".rds")))[[ref_period]]$ref_means %>%
  dplyr::select(pop,longitude,latitude,any_of(clim_var))

# Loading future climatic data extracted from the rasters for each GCM
list_clim_fut <- readRDS(here("data/ClimaticData/MaritimePinePops/ClimatePopulationLocationValuesExtractedFromRasters_FiveGCMs_2041-2070_SSP370.rds")) %>% 
  lapply(function(clim_fut){clim_fut %>% dplyr::select(pop,gcm,any_of(clim_var))})

Warning! (from the GDM website) Note that if your site coordinates are longitude-latitude, the calculation of geographic distances between sites will have errors, the size of which will depend on the geographic extent and location of your study region. We hope to deal with this in a later release, but for now you can avoid these problems by using a projected coordinate system (e.g., equidistant).

Therefore, we reproject the geographic (spheroid) CRS of the population coordinates (with units in degrees longitude and latitude) to a projected (two-dimensional; cartesian) CRS (typically with units of meters from a datum). We will use the CRS EPSG:3035.

Code
# Steps:
  # extract the geographic coordinates of the populations
  # transform the coordinates in spatial points and specify the CRS (WGS84) with sp package
  # transform to a SpatVector object for the terra package
  # reproject with terra package in CRS EPSG:3035
  # extract population coordinates from the SpatVector with terra package

pop_coord_proj <- clim_past %>%
  dplyr::select(contains("ude")) %>% 
  sp::SpatialPoints(proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")) %>% 
  terra::vect() %>% 
  terra::project("EPSG:3035") %>% 
  terra::crds() %>% 
  as_tibble() %>% 
  dplyr::rename(long_EPSG_3035=x, lat_EPSG_3035=y)

# merge the population coordinates with the climatic variables
clim_past <- bind_cols(clim_past,pop_coord_proj) %>% dplyr::select(-contains("ude"))
list_clim_fut <- lapply(list_clim_fut, function(clim_fut) bind_cols(clim_fut,pop_coord_proj) %>% dplyr::select(-contains("ude")))

3 GDM fitting

Combining genomic and climatic data with formatsitepair

We use the formatsitepair function to combine the genomic and climatic data into population-pair table format. In our case, we use a pairwise biological dissimilarity matrix (i.e. pairwise \(F_{ST}\)) as response variable, and therefore we have to set bioFormat=3 in the formatsitepair function.

Warning! The rows and columns of the distance matrix have to be in the same order as the rows of the climatic variables. And the distance matrix and the climatic data must not include NAs.

GDM fitting

From the GDM website: ’GDM is a nonlinear extension of permutational matrix regression that uses flexible splines and generalized linear modeling (GLM) to accommodate two types of nonlinearity common in ecological datasets: (1) variation in the rate of compositional turnover (non-stationarity) along environmental gradients, and (2) the curvilinear relationship between biological distance and environmental and geographical distance.

Generalized dissimilarity models are fitted with the gdm function.

Different options are available:

  • geo=T to specify that the model should be fit with geographical distance.

  • splines an optional vector specifying the the number of I-spline basis functions (the default is three, with larger values producing more complex splines).

  • knots an optional vector specifying the locations of “knots” along the splines (defaults to 0 (minimum), 50 (median), and 100 (maximum) quantiles when three I-spline basis functions are used). E

From the GDM website: ‘Even though these option are available, using the default values for these parameters will work fine for most applications. In other words, unless you have a good reason, you should probably use the default settings for splines and knots. The effects (and significance) of altering the number of splines and knot locations has not been systematically explored.’

Code
snp_sets <- lapply(snp_sets, function(x){
  
x$gdm_tab <- gdm_tab <- formatsitepair(bioData=fst_matrices[[x$set_code]], 
                                       bioFormat=3,
                                       XColumn="long_EPSG_3035", 
                                       YColumn="lat_EPSG_3035", 
                                       predData=clim_past, 
                                       siteColumn="pop")  

x$gdm_mod <- gdm_mod  <- gdm(data=gdm_tab, geo=TRUE)

return(x)

})

4 GDM plots

Observed genetic distance vs predicted climatic distance

Code
# Plot predicted climatic distance vs observed genetic distance

plots <- lapply(snp_sets, function(x){
  
overlay_x <- seq( from=min(x$gdm_mod$ecological), to=max(x$gdm_mod$ecological), length=length(x$gdm_mod$ecological) )
overlay_y <- 1 - exp( - overlay_x )
dfline <- tibble(x=overlay_x,y=overlay_y)  

tibble(clim=x$gdm_mod$ecological,
       obs=x$gdm_mod$observed) %>% 
  ggplot(aes(x=clim,y=obs)) + 
  geom_point(color=rgb(0,0,1,0.5)) +
  geom_line(data=dfline,aes(x=x,y=y), color="gray60") +
  xlim(c(0,1.6)) + 
  ylim(c(0,1)) + 
  xlab("Predicted climatic distance") +
  ylab("Observed genetic distance") +
  ggtitle(x$set_name) + 
  theme_bw() +
  theme(title = element_text(size=8))
  
})

big_plot <- plot_grid(plotlist = plots, nrow=2)

ggsave(filename = here("figs/GDM/PredictedClimaticDistancevsObservedGeneticDistancePlots.pdf"), 
       big_plot,
       device="pdf", width=10, height=6)

big_plot

Predicted versus observed genetic distance

Code
plots <- lapply(snp_sets, function(x){
  
tibble(obs=x$gdm_mod$observed, # can also be extracted with x$gdm_tab$distance
       pred=x$gdm_mod$predicted) %>% 
  ggplot(aes(x=obs,y=pred)) + 
  geom_abline(intercept = 0, slope = 1, color="gray60") +
  geom_point(color=rgb(0,0,1,0.5)) +
  xlim(c(0,1)) + ylim(c(0,1)) + 
  xlab("Observed genetic distance") +
  ylab("Predicted genetic distance") +
  ggtitle(x$set_name) + 
  theme_bw()+
  theme(title = element_text(size=8))
  
})

big_plot <- plot_grid(plotlist = plots, nrow=2)

ggsave(filename = here("figs/GDM/PredictedvsObservedGeneticDistancePlots.pdf"), 
       big_plot,
       device="pdf", width=10, height=6)

big_plot

I-splines

From the GDM website: ‘The fitted I-splines provide an indication of how population genetic composition changes along each climatic gradient. They are one of the most informative components of a fitted GDM and so plotting and scrutinizing the splines is a major part of interpreting GDM and the analyzed biological patterns.’

I-splines are shown for predictors with non-zero coefficients, i.e. predictors with a relationship with the genetic distance.

I-spline interpretation: (from the GDM website) ’The maximum height of each spline indicates the magnitude of total genetic change along that gradient and thereby corresponds to the relative importance of that predictor in contributing to allelic turnover while holding all other variables constant (i.e., is a partial climatic distance). The spline’s shape indicates how the rate of genetic change varies with position along that gradient. Thus, the splines provide insight into the total magnitude of genetic change as a function of each gradient and where along each gradient those changes are most pronounced.

Code
lapply(snp_sets, function(x){
  
spline_data <- isplineExtract(x$gdm_mod)

# extract variables with importance different from 0
predictors <- which(apply(spline_data$y,2, sum) != 0) %>% names() 

plots <- lapply(predictors, function(predictor){
 
if(predictor=="Geographic"){
  predictor_name <- "Geography (km)"
  spline_data$x[,predictor] <- spline_data$x[,predictor] / 1000
}  else {
predictor_name <- extract_climatedt_metadata(var_clim = clim_var) %>% 
  filter(label %in% predictor) %>% 
  mutate(predictor_legend= paste0(description, " (", label,"; ",unit_symbol,")")) %>% 
  pull(predictor_legend)
}
  
tibble(x=spline_data$x[,predictor],
       y=spline_data$y[,predictor]) %>% 
  ggplot() + 
  geom_line(aes(x=x,y=y),color="blue",linewidth=2) +
  ylim(c(0,1.3)) + 
  ylab("Partial genetic distance") +
  xlab(predictor_name) +
  theme(axis.title.x = element_text(size=9)) +
  theme_bw() 
  
})

# make title
title <- ggdraw() + 
  draw_label(x$set_name,
             fontface = 'bold',
             x = 0,
             hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

# merge plots
p <- plot_grid(plotlist=plots)

# If we want to save the figures without title
# ggsave(here(paste0("figs/GDM/IsplinePlot_",x$set_code,"_noTitle.pdf")), p, width=13,height=9, device="pdf")


# merge plots + title
p <- plot_grid(title, p, ncol = 1,rel_heights = c(0.1, 1))

# save the plots
ggsave(here(paste0("figs/GDM/IsplinePlot_",x$set_code,".pdf")), p, width=13,height=9, device="pdf")

})

5 GDM projections

Warning! Note that the formatsitepair function assumes that the coordinates of the sites are in the same coordinate system as the rasters. At present, no checking is performed to ensure this is the case..

6 GDM predictions

From the GDM website: A GDM model can be used to (i) predict the genetic dissimilarity between population pairs in space or between times using the predict function and (ii) transform the predictor variables from their arbitrary climatic scales to a common biological importance scale using the gdm.transform function.

6.1 Genomic offset of the populations

We want to predict the expected degree of maladaptation in units of the response variable (\(F_{ST}\)) for each population.

6.1.1 At the location of the populations

We aim to predict the genomic offset at the location of the studied populations (predictions based on spatial points). For that, we build a dataset with:

  • the genetic distance fixed to O. Indeed, we want the disruption of the current gene-climate relationships under future climates, so we assume that the genetic composition is the same between current and future climate to do this calculation.

  • weights are fixed to 1.

Code
# dataframe with past climatic data
pop_tab_pred_clim_past <- clim_past %>% 
  dplyr::select(-pop) %>% 
  dplyr::rename(xCoord=long_EPSG_3035,yCoord=lat_EPSG_3035) %>% 
  set_colnames(paste0("s1.",colnames(.)))

# list of dataframe with future climatic data for each GCM
list_pop_tab_pred <- list_clim_fut %>% 
  lapply(function(clim_fut){
  
clim_fut <-  clim_fut %>% 
  dplyr::select(-pop,-gcm) %>% 
  dplyr::rename(xCoord=long_EPSG_3035,yCoord=lat_EPSG_3035) %>% 
  set_colnames(paste0("s2.",colnames(.)))

bind_cols(pop_tab_pred_clim_past, clim_fut) %>% 
  mutate(distance=0,weights=0) %>% 
  dplyr::select(distance,weights,contains("Coord"), everything())})
  
# calculate the genomic offset for each GCM and each snp set
snp_sets <- lapply(snp_sets, function(x){

x$go <- lapply(list_pop_tab_pred, function(pop_tab_pred){
  
predict(x$gdm_mod,pop_tab_pred)

})
return(x)
})
Code
# Function to plot the relationship between the Euclidean climatic distance and the GDM genetic offset
source(here("scripts/functions/make_eucli_plot.R"))

# Calculate the Euclidean climatic distance
source(here("scripts/functions/generate_scaled_clim_datasets.R"))
clim_var <- readRDS(here("data/ClimaticData/NamesSelectedVariables.rds"))
clim_dfs <- generate_scaled_clim_datasets(clim_var,clim_ref_adj = FALSE)


# Calculate the Euclidean climatic distance
list_dist_env <- clim_dfs$clim_pred %>% lapply(function(clim_pred){
  
Delta = clim_dfs$clim_ref %>% dplyr::select(any_of(clim_var)) - clim_pred %>% dplyr::select(any_of(clim_var)) 
dist_env = sqrt( rowSums(Delta^2) )

})

# list_dist_env <- list_clim_fut %>% lapply(function(clim_fut){
#   
# Delta = clim_past %>% dplyr::select(-pop,-contains("EPSG")) - clim_fut %>% dplyr::select(-pop,-gcm,-contains("EPSG")) 
# dist_env = sqrt( rowSums(Delta^2) )
# 
# })

# Main gene pools (for the figures)
gps <- readRDS(here("data/GenomicData/MainGenePoolPopulations.rds")) %>% 
  arrange(pop)
Code
par(mfrow=c(1,2))


lapply(snp_sets, function(x) {

lapply(names(list_dist_env), function(gcm){
  
make_eucli_plot(
  X = list_dist_env[[gcm]],
  Y = x$go[[gcm]],
  colors = gps$color_main_gp_pop,
  color_names = gps$main_gp_pop,
  ylab = "GDM genomic offset",
  legend_position="topright",
  plot_title = paste0(x$set_name," - ", gcm))

})
}) 

Code
# We generate scatter plots for the Supplementary Information.
# ============================================================

# Axis limits
# ===========
max_go <- lapply(snp_sets, function(z){
  z$go %>% unlist()
}) %>% unlist() %>% max()

range_eucli <- list_dist_env %>% unlist() %>% range()

# Run the function
# ================
lapply(snp_sets, function(set_i) {

p <- lapply(names(list_dist_env), function(gcm){
  
make_ggscatterplot(
  x = list_dist_env[[gcm]],
  y = set_i$go[[gcm]],
  title = gcm,
  range_eucli = range_eucli,
  max_go = max_go)

})

# remove y-labels to graphs in the second column
p[[2]] <- p[[2]] + ylab("")
p[[4]] <- p[[4]] + ylab("")


# remove x-labels to graphs in the second and third rows
p[[1]] <- p[[1]] + xlab("")
p[[2]] <- p[[2]] + xlab("")
p[[3]] <- p[[3]] + xlab("")

p[[6]] <- get_legend(p[[1]])

for(i in 1:5){p[[i]] <- p[[i]]  +  theme(legend.position = "none")} 


p_grid <- plot_grid(plotlist=p, nrow = 3) 

# We save the figure in pdf
p_grid %>% 
  ggsave(here(paste0("figs/GDM/ScatterPlotEucliDistance_",set_i$set_code,".pdf")), 
         .,
         width=7,
         height=8,
         device="pdf")

# we save the figure in png
p_grid %>% 
  ggsave(here(paste0("figs/GDM/ScatterPlotEucliDistance_",set_i$set_code,".png")), 
         .,
         width=7,
         height=8)

})
$all_cand
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/GDM/ScatterPlotEucliDistance_all_cand.png"

$common_cand
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/GDM/ScatterPlotEucliDistance_common_cand.png"

$corrected_cand
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/GDM/ScatterPlotEucliDistance_corrected_cand.png"

$randomfreq_control
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/GDM/ScatterPlotEucliDistance_randomfreq_control.png"

$simfreq_control
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/GDM/ScatterPlotEucliDistance_simfreq_control.png"

$all
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/GDM/ScatterPlotEucliDistance_all.png"
Code
# Function to make the genomic offset maps
source(here("scripts/functions/make_go_map.R"))

# Population coordinates
pop_coord <-  readRDS(here(paste0("data/ClimaticData/MaritimePinePops/ClimatePopulationLocationPointEstimates_ReferencePeriods_",adj,".rds")))[[ref_period]]$ref_means %>%
  dplyr::select(pop,longitude,latitude)


# Generate the maps for each set of SNPs and each GCM
lapply(snp_sets, function(x) {

go_maps <- lapply(names(list_clim_fut), function(gcm){
  
df <- pop_coord %>% mutate(GO = x$go[[gcm]])

make_go_map(df = df,
            plot_title = gcm,
            go_limits = c(0, max(x$go[[gcm]])),
            legend_box_background = "white",
            point_size = 3)

})

legend_maps  <- get_legend(go_maps[[1]])

go_maps <- lapply(go_maps, function(y) y + theme(legend.position = "none"))

go_maps$legend_maps <- legend_maps

go_maps <-plot_grid(plotlist=go_maps)

ggsave(here(paste0("figs/GDM/GOMaps_PopLocations_",x$set_code,".pdf")), go_maps, width=10,height=6, device="pdf")
ggsave(here(paste0("figs/GDM/GOMaps_PopLocations_",x$set_code,".png")), go_maps, width=10,height=6)


# =========
# Add title
# =========
title <- ggdraw() + 
  draw_label(
    x$set_name,
    fontface = 'bold',
    x = 0,
    hjust = 0
  ) +
  theme(plot.margin = margin(0, 0, 0, 7))

# merge title and plots
plot_grid(
  title, go_maps,
  ncol = 1,
  rel_heights = c(0.1, 1))

  })

For each GCM, we attribute the value 1 to the top five populations with the highest genomic offset and we attribute the value 0 to the other populations. We then count the number of 1 for each population, which gives the table and map below:

Code
selected_SNP_set <- "common_cand"

source(here("scripts/functions/make_high_go_pop_maps.R"))

high_go_pops <- make_high_go_pop_maps(pop_coord=pop_coord,
                                      list_go = snp_sets[[selected_SNP_set]]$go,
                                      ggtitle = paste0("GDM - ",snp_sets[[selected_SNP_set]]$set_name),
                                      nb_id_pop = 5) # number of selected populations
Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2 3.5.0.
ℹ Please use the `legend.position.inside` argument of `theme()` instead.
Code
saveRDS(high_go_pops, file = here("outputs/GDM/high_go_pops.rds"))

high_go_pops[[1]] %>% kable_mydf
pop longitude latitude GFDL-ESM4 IPSL-CM6A-LR MPI-ESM1-2-HR MRI-ESM2-0 UKESM1-0-LL sum_go
VAL -4.31 40.52 0 0 0 0 0 0
ARN -5.12 40.19 0 0 0 0 0 0
CEN -4.49 40.28 0 0 0 0 0 0
BON -1.66 39.99 0 0 0 0 0 0
ORI -2.35 37.53 0 0 0 0 0 0
QUA -0.36 38.97 0 0 0 0 0 0
OLB -0.62 40.17 0 0 0 0 0 0
PIA 9.46 42.02 0 0 0 0 0 0
VER -1.09 45.55 0 0 0 0 0 0
MAD -5.23 35.18 0 0 0 0 0 0
COC -4.50 41.25 0 0 0 0 0 0
CUE -4.48 41.37 0 0 0 0 0 0
COM -3.95 36.83 0 0 0 0 0 0
TAM -5.02 33.60 0 0 0 0 0 0
HOU -1.15 45.18 0 0 0 0 0 0
BAY -2.88 41.52 0 0 0 0 0 0
CAR -4.28 41.17 0 0 0 0 0 0
STJ -2.03 46.76 0 0 0 0 0 0
PIE 9.04 41.97 0 0 0 0 0 0
OLO -1.83 46.57 0 0 0 0 0 0
MIM -1.30 44.13 0 0 0 0 0 0
PLE -2.34 47.78 0 0 0 0 0 0
SAL -3.06 41.83 0 0 0 0 0 0
SIE -6.49 43.53 0 0 0 0 1 1
LEI -8.96 39.78 0 0 0 1 0 1
PET -1.30 44.06 0 1 0 0 0 1
CAS -6.98 43.50 0 0 0 0 1 1
PUE -6.63 43.55 0 0 0 0 1 1
LAM -6.22 43.56 0 1 1 0 0 2
SEG -8.45 42.82 1 0 0 1 0 2
SAC -8.36 42.12 1 0 1 1 0 3
CAD -6.42 43.54 1 1 1 0 1 4
ALT -6.49 43.28 1 1 1 1 0 4
ARM -6.46 43.30 1 1 1 1 1 5
Code
high_go_pops[[2]]

6.1.1.1 Comparing GO predictions

We look at the correlation across the different genomic offset predictions at the location of the populations, i.e. those based on all SNPs and those based on sets of candidates or control SNPs.

Code
lapply(names(snp_sets[[1]]$go), function(gcm){
  
lapply(snp_sets, function(x) x$go[[gcm]]) %>% 
  as_tibble() %>%
  cor() %>% 
  corrplot(method = 'number',type = 'lower', 
           diag = FALSE,mar=c(0,0,2,0),
           title=gcm,
           number.cex=2,tl.cex=1.5)
  
})

6.1.2 Across the species range

To predict the genomic offset across the species range, we use climatic rasters for the reference period (1901-1950). These rasters were provided by Maurizio Marchi and come from the Climate Downscaling Tool (ClimateDT). The values of climatic variables at the location of the populations match between the rasters and the spatial points (I checked it in the report 3_CheckingPastFutureClimatesPopulationLocations).

6.1.2.1 Maps with raster package

Code
# Past rasters
path <- here("data/ClimaticData/ClimateDTRasters/1km_1901-1950_Extent-JulietteA/")
tif_paths <- lapply(clim_var, function(x) paste0(path,"/",x,".tif"))
past_rasts <- raster::stack(tif_paths) %>% projectRaster(crs = "EPSG:3035")


list_pred_rasts <- lapply(snp_sets, function(x){ # for each snp set
lapply(names(list_clim_fut), function(gcm){ # for each GCM
  
# Future rasters  
path <- here(paste0("data/ClimaticData/ClimateDTRasters/1km_",gcm,"_2041-2070_ssp370_Extent-JulietteA/"))
tif_paths <- lapply(clim_var, function(x) paste0(path,"/",x,".tif"))
fut_rasts <- raster::stack(tif_paths)
fut_rasts <- projectRaster(fut_rasts, crs = "EPSG:3035")
  

predict(x$gdm_mod,
        past_rasts, # raster with current climate at 2.5 minutes resolution
        time=TRUE,
        fut_rasts)  # Rasters with future climates at 2.5 minutes resolution

}) %>% setNames(names(list_clim_fut))
})


list_pred_rasts %>% saveRDS(here("outputs/GDM/go_pred_rasters.rds"))
Code
list_pred_rasts <- readRDS(here("outputs/GDM/go_pred_rasters.rds"))
Code
par(mfrow=c(1,2))

names(snp_sets) %>% lapply(function(x){
  
lapply(names(list_clim_fut), function(gcm){

raster::plot(list_pred_rasts[[x]][[gcm]], 
             col=rgb.tables(1000),
             axes=FALSE, box=FALSE,
             main=paste0(snp_sets[[x]]$set_name," - ", gcm))
})
  
})

6.1.2.2 Maps with ggplot2

We project the genomic offset predictions for the set of candidate SNPs and the mean genomic offset across the five GCMs.

MANUSCRIPT FIGURE: This figure corresponds to Figure 6a in the main manuscript.

Code
# Buffer for the maps (maritime pine distribution)
range_buffer = shapefile(here('data/Mapping/PinpinDistriEUforgen_NFIplotsBuffer10km.shp'))

# We project the genomic offset only for the candidate SNPs
snp_set <- snp_sets[["all_cand"]]

# We load the rasters with the climates of the reference period
path <- here("data/ClimaticData/ClimateDTRasters/1km_1901-1950_Extent-JulietteA/")
ref_rasts <- lapply(clim_var, function(x) paste0(path,"/",x,".tif")) %>% 
  raster::stack() %>%
  mask(range_buffer) %>% 
  projectRaster(crs = "EPSG:3035")


df <- lapply(names(list_clim_fut), function(gcm){ # for each GCM
  
# Future rasters  
path <- here(paste0("data/ClimaticData/ClimateDTRasters/1km_",gcm,"_2041-2070_ssp370_Extent-JulietteA/"))
fut_rasts <- lapply(clim_var, function(x) paste0(path,"/",x,".tif")) %>% 
  raster::stack() %>% 
  mask(range_buffer) %>% 
  projectRaster(crs = "EPSG:3035")
  

predict(snp_set$gdm_mod,
              ref_rasts, # rasters with reference climate at 2.5 minutes resolution
              time=TRUE,
              fut_rasts) %>%  # rasters with future climates at 2.5 minutes resolutio
  projectRaster(crs = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0") %>% 
  rasterToPoints() %>% 
  set_colnames(c("x","y","GO")) %>% 
  as_tibble() %>% 
  mutate(gcm=gcm)

}) %>% 
  bind_rows %>% 
  pivot_wider(values_from = GO, names_from = gcm) %>% 
  mutate(mean_GO = rowMeans(dplyr::select(.,-c(x:y))))


df %>% saveRDS(here("outputs/GDM/go_pred_rasters_2_allcand.rds"))
Code
# Map options
# ===========
point_size = 2
x_limits = c(-10, 15)
y_limits = c(31, 52)

# Country borders
world <- ne_countries(scale = "medium", returnclass = "sf")

# Load the mean GO projections
df <- readRDS(here("outputs/GDM/go_pred_rasters_2_allcand.rds"))

p <- ggplot(data=df) + 
  geom_sf(data = world, fill="gray98") + 
  scale_x_continuous(limits = x_limits) +
  scale_y_continuous(limits = y_limits) + 
  geom_raster(aes(x = x, y = y, fill = mean_GO), alpha = 1) + 
  scale_fill_gradient2(low="blue", mid= "yellow", high="red",
                       midpoint=(max(df$mean_GO)-min(df$mean_GO))/2,
                       limits=c(min(df$mean_GO),max(df$mean_GO)),
                       name = "Genomic offset from GDM") +
  xlab("Longitude") + ylab("Latitude") +
  #ggtitle("Genomic offset predictions with GDM") +
  theme_bw(base_size = 11) +
  theme(panel.grid = element_blank(), 
        plot.background = element_blank(), 
        panel.background = element_blank(), 
        legend.position = c(0.8,0.15),
        legend.box.background = element_rect(colour = "gray80"),
        legend.title = element_text(size=10),
        strip.text = element_text(size=11))

p %>% ggsave(here("figs/GDM/GOmeanProjections_AllCandSNPs.pdf"),., width=6,height=6, device="pdf")
p %>% ggsave(here("figs/GDM/GOmeanProjections_AllCandSNPs.png"),., width=6,height=6)
p

6.1.3 Corr btw predictions with rasters or spatial points

We check that the genomic offset predictions at the locations of the populations are correlated between those obtained with the spatial points and those obtained with the rasters.

Code
# checking correlations
names(snp_sets) %>% lapply(function(x){
  
cor_go <- lapply(names(list_clim_fut), function(gcm){
  
go_rast <- raster::extract(list_pred_rasts[[x]][[gcm]], clim_past[,c("long_EPSG_3035","lat_EPSG_3035")])

cor(snp_sets[[x]]$go[[gcm]],go_rast) 
  
}) %>% unlist() 
  
tibble(gcm=names(list_clim_fut),
            cor_go=cor_go)
  
}) %>% 
  setNames(names(snp_sets)) %>% 
  bind_rows(.id="snp_set") %>% 
  kable_mydf()
snp_set gcm cor_go
all_cand GFDL-ESM4 0.99
all_cand IPSL-CM6A-LR 0.99
all_cand MPI-ESM1-2-HR 0.99
all_cand MRI-ESM2-0 0.99
all_cand UKESM1-0-LL 0.99
common_cand GFDL-ESM4 0.99
common_cand IPSL-CM6A-LR 0.99
common_cand MPI-ESM1-2-HR 0.99
common_cand MRI-ESM2-0 0.99
common_cand UKESM1-0-LL 0.99
corrected_cand GFDL-ESM4 0.99
corrected_cand IPSL-CM6A-LR 0.99
corrected_cand MPI-ESM1-2-HR 0.99
corrected_cand MRI-ESM2-0 0.99
corrected_cand UKESM1-0-LL 0.99
randomfreq_control GFDL-ESM4 1.00
randomfreq_control IPSL-CM6A-LR 1.00
randomfreq_control MPI-ESM1-2-HR 1.00
randomfreq_control MRI-ESM2-0 1.00
randomfreq_control UKESM1-0-LL 1.00
simfreq_control GFDL-ESM4 0.99
simfreq_control IPSL-CM6A-LR 0.99
simfreq_control MPI-ESM1-2-HR 0.99
simfreq_control MRI-ESM2-0 1.00
simfreq_control UKESM1-0-LL 0.99
all GFDL-ESM4 0.99
all IPSL-CM6A-LR 0.99
all MPI-ESM1-2-HR 1.00
all MRI-ESM2-0 1.00
all UKESM1-0-LL 1.00

This is ok!

7 Validation - NFI plots

Climatic data and coordinates

Code
# Load the climatic data of the NFI plots.
nfi_clim <- readRDS(here("data/ClimaticData/NFIplots/NFIclimate.rds")) 

# Reproject the plotcode coordinates in the projection system `EPSG:3035`.
plotcoord_EPSG_3035 <- nfi_clim[[1]] %>%
  dplyr::select(contains("ude")) %>% 
  sp::SpatialPoints(proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")) %>% 
  terra::vect() %>% 
  terra::project("EPSG:3035") %>% 
  terra::crds() %>% 
  as_tibble() %>% 
  dplyr::rename(xCoord=x, yCoord=y)

# Keep only the climatic variables of interest and the plotcode coordinates
nfi_clim <- nfi_clim %>% 
  lapply(function(x){
    
x %>% bind_cols(plotcoord_EPSG_3035) %>% 
      dplyr::select(contains("Coord"),any_of(clim_var))})

# Format as GDM inputs
nfi_clim[[1]] <- nfi_clim[[1]] %>% set_colnames(paste0("s1.",colnames(.)))
nfi_clim[[2]] <- nfi_clim[[2]] %>% set_colnames(paste0("s2.",colnames(.)))

nfi_clim <- bind_cols(nfi_clim) %>% 
  mutate(distance=0,weights=0) %>% 
  dplyr::select(distance,weights,contains("Coord"), everything())

Calculate the genomic offset for each snp set

Code
snp_sets <- lapply(snp_sets, function(x){

x$go_nfi <- predict(x$gdm_mod,nfi_clim)

return(x)
})

Map genomic offset predictions at the location of the NFI plots

These figures are included in the Supplementary Information.

Code
source(here("scripts/functions/make_go_map.R"))

lapply(snp_sets, function(x) {

df <- readRDS(here("data/ClimaticData/NFIplots/NFIclimate.rds"))[[1]] %>% 
dplyr::select(contains("ude")) %>% mutate(GO = x$go_nfi) 

p <- make_go_map(df = df, 
                 type="NFI",
                 point_size = 0.5,
                 plot_title =  x$set_name,
                 go_limits = c(0,max(x$go_nfi)),
                 legend_position = c(0.85,0.2), 
                 y_limits = c(35, 51))
  
  
ggsave(here(paste0("figs/GDM/NFI_GOmap_",x$set_code,".pdf")), p, width=6,height=6, device="pdf")
ggsave(here(paste0("figs/GDM/NFI_GOmap_",x$set_code,".png")), p, width=6,height=6)
# p <- p + theme(plot.title = element_blank()) # to remove the title

# Show maps in the Quarto document
# ================================
p  
  })

We look at the correlation across the different genomic offset predictions in the NFI plots, i.e. those based on all SNPs and those based on sets of candidates or control SNPs.

Code
lapply(snp_sets, function(x) x$go_nfi) %>% 
  as_tibble() %>%
  cor() %>% 
  corrplot(method = 'number',type = 'lower', diag = FALSE,mar=c(0,0,2,0),
               number.cex=2,tl.cex=1.5)

8 Validation - Common gardens

Climatic data and coordinates

Code
# format climatic data
cg_clim <- readRDS(here("data/ClimaticData/CommonGardens/ClimateCG.rds")) %>%  
  dplyr::select(any_of(clim_var)) %>% 
  set_colnames(paste0("s2.",colnames(.))) %>%  
  replicate(nrow(pop_tab_pred_clim_past),.,simplify=F) %>% 
  bind_rows() %>% 
  bind_cols(slice(pop_tab_pred_clim_past, rep(1:n(), each = 5))) %>% 
  mutate(s2.xCoord=s1.xCoord, # same coordinates, so that only climatic differences are considered
         s2.yCoord=s1.yCoord,
         distance=0,
         weights=0) %>% 
  dplyr::select(distance, weights, s1.xCoord, s1.yCoord, s2.xCoord, s2.yCoord,contains("s1"), contains("s2"))

Calculate the genomic offset for each snp set

Code
go_df <- readRDS(here("data/ClimaticData/CommonGardens/ClimateCG.rds")) %>% 
  dplyr::select(cg) %>% 
  replicate(nrow(pop_tab_pred_clim_past),.,simplify=F) %>% 
  bind_rows() %>% 
  bind_cols(slice(clim_past[,"pop"],rep(1:n(), each = 5)))

snp_sets <- lapply(snp_sets, function(x){

x$go_cg <- go_df %>% 
  mutate(go=predict(x$gdm_mod,cg_clim)) %>% 
  pivot_wider(names_from = cg, values_from = go)

return(x)
})

Map genomic offset predictions at the locations of the populations

Code
cg_coord <- readRDS(here("data/ClimaticData/CommonGardens/ClimateCG.rds")) %>% dplyr::select(cg,contains("ude"))
cg_names <- unique(cg_coord$cg)

go_maps_cg <- lapply(cg_names, function(cg_name){

p <- lapply(snp_sets, function(x) {

df <- pop_coord %>%
        left_join(x$go_cg[,c("pop",cg_name)], by="pop") %>% 
        dplyr::rename(GO=all_of(cg_name))
  
 p <- make_go_map(df = df,
                  plot_title = paste0(str_to_title(cg_name), " - ",x$set_name),
                  point_size = 3,
                  type="CG",
                  go_limits = c(0,max(x$go_cg[[cg_name]])),
                  cg_coord = filter(cg_coord, cg == cg_name))
 
 
ggsave(filename = here(paste0("figs/GDM/GOmap_",x$set_code,"_",cg_name,".pdf")), device = "pdf",width=6,height=6)
ggsave(filename = here(paste0("figs/GDM/GOmap_",x$set_code,"_",cg_name,".png")), width=6,height=6)

# p + theme(plot.title = element_blank(), legend.position = "none") # to rm the title
 
p
 
  })

plot_grid(plotlist=p,nrow=2)
  
}) %>% setNames(cg_names)

pdf(here("figs/GDM/GOmaps_CGs.pdf"), width=15,height=6)
lapply(go_maps_cg, function(x) x)
dev.off()

# show maps
lapply(go_maps_cg, function(x) x)

We look at the correlation across the different genomic offset predictions in the common gardens, i.e. those based on all SNPs and those based on sets of candidates or control SNPs.

Code
lapply(cg_names, function(cg_name){

lapply(snp_sets, function(x) x$go_cg) %>% 

lapply(function(set){
  set[[cg_name]]
}) %>% 
    as_tibble() %>% 
    cor() %>% 
    corrplot(method = 'number',type = 'lower', 
             diag = FALSE,mar=c(0,0,2,0),
             title=str_to_title(cg_name),
             number.cex=2,tl.cex=1.5)

})

9 Saving GO predictions

Let’s save the genomic offset predictions for comparison with the other methods.

Code
snp_sets %>% saveRDS(file=here("outputs/GDM/go_predictions.rds"))

10 Session information

Code
devtools::session_info()
─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.4.2 (2024-10-31)
 os       Ubuntu 22.04.5 LTS
 system   x86_64, linux-gnu
 ui       X11
 language (EN)
 collate  fr_FR.UTF-8
 ctype    fr_FR.UTF-8
 tz       Europe/Paris
 date     2025-01-16
 pandoc   3.1.11 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/x86_64/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
 package           * version   date (UTC) lib source
 cachem              1.0.8     2023-05-01 [1] CRAN (R 4.4.0)
 cellranger          1.1.0     2016-07-27 [1] CRAN (R 4.4.0)
 class               7.3-23    2025-01-01 [4] CRAN (R 4.4.2)
 classInt            0.4-10    2023-09-05 [1] CRAN (R 4.4.0)
 cli                 3.6.2     2023-12-11 [1] CRAN (R 4.4.0)
 cluster             2.1.8     2024-12-11 [4] CRAN (R 4.4.2)
 codetools           0.2-19    2023-02-01 [4] CRAN (R 4.2.2)
 colorRamps        * 2.3.4     2024-03-07 [1] CRAN (R 4.4.0)
 colorspace          2.1-0     2023-01-23 [1] CRAN (R 4.4.0)
 corrplot          * 0.92      2021-11-18 [1] CRAN (R 4.4.0)
 cowplot           * 1.1.3     2024-01-22 [1] CRAN (R 4.4.0)
 DBI                 1.2.2     2024-02-16 [1] CRAN (R 4.4.0)
 devtools            2.4.5     2022-10-11 [1] CRAN (R 4.4.0)
 digest              0.6.35    2024-03-11 [1] CRAN (R 4.4.0)
 doParallel          1.0.17    2022-02-07 [1] CRAN (R 4.4.0)
 dplyr             * 1.1.4     2023-11-17 [1] CRAN (R 4.4.0)
 e1071               1.7-14    2023-12-06 [1] CRAN (R 4.4.0)
 ellipsis            0.3.2     2021-04-29 [1] CRAN (R 4.4.0)
 evaluate            0.23      2023-11-01 [1] CRAN (R 4.4.0)
 fansi               1.0.6     2023-12-08 [1] CRAN (R 4.4.0)
 farver              2.1.2     2024-05-13 [1] CRAN (R 4.4.0)
 fastmap             1.1.1     2023-02-24 [1] CRAN (R 4.4.0)
 forcats           * 1.0.0     2023-01-29 [1] CRAN (R 4.4.0)
 foreach             1.5.2     2022-02-02 [1] CRAN (R 4.4.0)
 fs                  1.6.4     2024-04-25 [1] CRAN (R 4.4.0)
 gdm               * 1.5.0-9.1 2022-12-01 [1] CRAN (R 4.4.0)
 generics            0.1.3     2022-07-05 [1] CRAN (R 4.4.0)
 ggplot2           * 3.5.1     2024-04-23 [1] CRAN (R 4.4.0)
 glue                1.7.0     2024-01-09 [1] CRAN (R 4.4.0)
 gtable              0.3.5     2024-04-22 [1] CRAN (R 4.4.0)
 here              * 1.0.1     2020-12-13 [1] CRAN (R 4.4.0)
 hierfstat         * 0.5-11    2022-05-05 [1] CRAN (R 4.4.0)
 highr               0.10      2022-12-22 [1] CRAN (R 4.4.0)
 hms                 1.1.3     2023-03-21 [1] CRAN (R 4.4.0)
 htmltools           0.5.8.1   2024-04-04 [1] CRAN (R 4.4.0)
 htmlwidgets         1.6.4     2023-12-06 [1] CRAN (R 4.4.0)
 httpuv              1.6.15    2024-03-26 [1] CRAN (R 4.4.0)
 httr                1.4.7     2023-08-15 [1] CRAN (R 4.4.0)
 iterators           1.0.14    2022-02-05 [1] CRAN (R 4.4.0)
 jsonlite            1.8.8     2023-12-04 [1] CRAN (R 4.4.0)
 kableExtra        * 1.4.0     2024-01-24 [1] CRAN (R 4.4.0)
 KernSmooth          2.23-26   2025-01-01 [4] CRAN (R 4.4.2)
 knitr             * 1.46      2024-04-06 [1] CRAN (R 4.4.0)
 labeling            0.4.3     2023-08-29 [1] CRAN (R 4.4.0)
 later               1.3.2     2023-12-06 [1] CRAN (R 4.4.0)
 latex2exp         * 0.9.6     2022-11-28 [1] CRAN (R 4.4.0)
 lattice             0.22-5    2023-10-24 [4] CRAN (R 4.3.1)
 lifecycle           1.0.4     2023-11-07 [1] CRAN (R 4.4.0)
 lubridate         * 1.9.3     2023-09-27 [1] CRAN (R 4.4.0)
 magrittr          * 2.0.3     2022-03-30 [1] CRAN (R 4.4.0)
 MASS                7.3-64    2025-01-04 [4] CRAN (R 4.4.2)
 Matrix              1.7-1     2024-10-18 [4] CRAN (R 4.4.1)
 memoise             2.0.1     2021-11-26 [1] CRAN (R 4.4.0)
 mgcv                1.9-1     2023-12-21 [4] CRAN (R 4.3.2)
 mime                0.12      2021-09-28 [1] CRAN (R 4.4.0)
 miniUI              0.1.1.1   2018-05-18 [1] CRAN (R 4.4.0)
 munsell             0.5.1     2024-04-01 [1] CRAN (R 4.4.0)
 nlme                3.1-166   2024-08-14 [4] CRAN (R 4.4.2)
 pbapply             1.7-2     2023-06-27 [1] CRAN (R 4.4.0)
 permute             0.9-7     2022-01-27 [1] CRAN (R 4.4.0)
 pillar              1.9.0     2023-03-22 [1] CRAN (R 4.4.0)
 pkgbuild            1.4.4     2024-03-17 [1] CRAN (R 4.4.0)
 pkgconfig           2.0.3     2019-09-22 [1] CRAN (R 4.4.0)
 pkgload             1.3.4     2024-01-16 [1] CRAN (R 4.4.0)
 plyr                1.8.9     2023-10-02 [1] CRAN (R 4.4.0)
 profvis             0.3.8     2023-05-02 [1] CRAN (R 4.4.0)
 promises            1.3.0     2024-04-05 [1] CRAN (R 4.4.0)
 proxy               0.4-27    2022-06-09 [1] CRAN (R 4.4.0)
 purrr             * 1.0.2     2023-08-10 [1] CRAN (R 4.4.0)
 R6                  2.5.1     2021-08-19 [1] CRAN (R 4.4.0)
 ragg                1.3.0     2024-03-13 [1] CRAN (R 4.4.0)
 raster            * 3.6-26    2023-10-14 [1] CRAN (R 4.4.0)
 Rcpp                1.0.12    2024-01-09 [1] CRAN (R 4.4.0)
 readr             * 2.1.5     2024-01-10 [1] CRAN (R 4.4.0)
 readxl            * 1.4.3     2023-07-06 [1] CRAN (R 4.4.0)
 remotes             2.5.0     2024-03-17 [1] CRAN (R 4.4.0)
 reshape2            1.4.4     2020-04-09 [1] CRAN (R 4.4.0)
 rlang               1.1.4     2024-06-04 [1] CRAN (R 4.4.0)
 rmarkdown           2.26      2024-03-05 [1] CRAN (R 4.4.0)
 rnaturalearth     * 1.0.1     2023-12-15 [1] CRAN (R 4.4.0)
 rnaturalearthdata * 1.0.0     2024-02-09 [1] CRAN (R 4.4.0)
 rprojroot           2.0.4     2023-11-05 [1] CRAN (R 4.4.0)
 rstudioapi          0.16.0    2024-03-24 [1] CRAN (R 4.4.0)
 scales              1.3.0     2023-11-28 [1] CRAN (R 4.4.0)
 sessioninfo         1.2.2     2021-12-06 [1] CRAN (R 4.4.0)
 sf                * 1.0-16    2024-03-24 [1] CRAN (R 4.4.0)
 shiny               1.8.1.1   2024-04-02 [1] CRAN (R 4.4.0)
 sp                * 2.1-4     2024-04-30 [1] CRAN (R 4.4.0)
 stringi             1.8.4     2024-05-06 [1] CRAN (R 4.4.0)
 stringr           * 1.5.1     2023-11-14 [1] CRAN (R 4.4.0)
 svglite             2.1.3     2023-12-08 [1] CRAN (R 4.4.0)
 systemfonts         1.0.6     2024-03-07 [1] CRAN (R 4.4.0)
 terra             * 1.7-71    2024-01-31 [1] CRAN (R 4.4.0)
 textshaping         0.3.7     2023-10-09 [1] CRAN (R 4.4.0)
 tibble            * 3.2.1     2023-03-20 [1] CRAN (R 4.4.0)
 tidyr             * 1.3.1     2024-01-24 [1] CRAN (R 4.4.0)
 tidyselect          1.2.1     2024-03-11 [1] CRAN (R 4.4.0)
 tidyverse         * 2.0.0     2023-02-22 [1] CRAN (R 4.4.0)
 timechange          0.3.0     2024-01-18 [1] CRAN (R 4.4.0)
 tzdb                0.4.0     2023-05-12 [1] CRAN (R 4.4.0)
 units               0.8-5     2023-11-28 [1] CRAN (R 4.4.0)
 urlchecker          1.0.1     2021-11-30 [1] CRAN (R 4.4.0)
 usethis             2.2.3     2024-02-19 [1] CRAN (R 4.4.0)
 utf8                1.2.4     2023-10-22 [1] CRAN (R 4.4.0)
 vctrs               0.6.5     2023-12-01 [1] CRAN (R 4.4.0)
 vegan               2.6-4     2022-10-11 [1] CRAN (R 4.4.0)
 viridisLite         0.4.2     2023-05-02 [1] CRAN (R 4.4.0)
 withr               3.0.0     2024-01-16 [1] CRAN (R 4.4.0)
 xfun                0.43      2024-03-25 [1] CRAN (R 4.4.0)
 xml2                1.3.6     2023-12-04 [1] CRAN (R 4.4.0)
 xtable            * 1.8-4     2019-04-21 [1] CRAN (R 4.4.0)
 yaml                2.3.8     2023-12-11 [1] CRAN (R 4.4.0)

 [1] /home/juliette/R/x86_64-pc-linux-gnu-library/4.4
 [2] /usr/local/lib/R/site-library
 [3] /usr/lib/R/site-library
 [4] /usr/lib/R/library

──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

References

Ferrier, Simon, Glenn Manion, Jane Elith, and Karen Richardson. 2007. “Using Generalized Dissimilarity Modelling to Analyse and Predict Patterns of Beta Diversity in Regional Biodiversity Assessment.” Diversity and Distributions 13 (3): 252–64. https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1472-4642.2007.00341.x.
Fitzpatrick, Matthew C, and Stephen R Keller. 2015. “Ecological Genomics Meets Community-Level Modelling of Biodiversity: Mapping the Genomic Landscape of Current and Future Environmental Adaptation.” Ecology Letters 18 (1): 1–16. https://onlinelibrary.wiley.com/doi/abs/10.1111/ele.12376.
Gougherty, Andrew V, Stephen R Keller, and Matthew C Fitzpatrick. 2021. “Maladaptation, Migration and Extirpation Fuel Climate Change Risk in a Forest Tree Species.” Nature Climate Change 11 (2): 166–71.
Mokany, Karel, Chris Ware, Skipton NC Woolley, Simon Ferrier, and Matthew C Fitzpatrick. 2022. “A Working Guide to Harnessing Generalized Dissimilarity Modelling for Biodiversity Analysis and Conservation Assessment.” Global Ecology and Biogeography 31 (4): 802–21. https://onlinelibrary.wiley.com/doi/full/10.1111/geb.13459.